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Abstract 


A numerical simulation tool for calculating the planar solid oxide fuel cells was described. The finite volume method was employed for 
the simulation, which was on the basis of the fundamental conservation laws of mass, momentum, energy and electrical charge. Temperature 
distributions, molar concentrations of gaseous species, current density and over potential were calculated using a single cell unit model with double 
channels of co-flow and counter-flow cases. The influences of operating conditions and anode structure on the performances of SOFC were also 
discussed. Simulation results show that the co-flow case has more uniform temperature and current density distributions and smaller temperature 
gradients, thus offers thermostructural advantages than the counter-flow case. Moreover, in co-flow case, with the increasing of delivery rate, 
temperature and hydrogen mass fraction of fuel, average temperature of PEN, current density and activation potential also rise. However, with 
increasing the delivery rate of air, average temperature of PEN decreases. In particular, it is effective to improve the output voltage by reducing the 


thickness of anode or increasing its porosity. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


The solid oxide fuel cell (SOFC) has been drawn attention 
because of its higher energy conversion efficiency, power den- 
sity, low environmental hazards and production cost. Thus, the 
SOFC is expected to reach commercialization in a few years and 
could be a promising alternative energy source for residential and 
distributed power plants in the 21st century [1-3]. However, fur- 
ther development of the planar SOFC faces challenges related to 
maximizing the power density and minimizing the non-uniform 
distribution of temperature, which contributes to thermal stress 
in the SOFC components [4-5]. 

The temperature and current density distributions in a SOFC 
are determined by the working conditions, such as the delivery 
rates and temperatures of the fuel and air to the system. Fuel 
utilization and fuel distribution are also critical due to the 
exothermic electrochemical reactions. Increased fuel flow tends 
to increase the uniformity of the reaction rates across the active 
area and decrease fuel utilization. Decreased fuel flow tends to 
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increase the fuel utilization, but can cause local fuel depletion 
and cold spots that exacerbate temperature non-uniformities. 
Air and fuel inlet temperatures also affect the reaction rates, 
cell temperature and fuel utilization. Therefore, management of 
the flow of air and fuel and the distribution of each, is critical to 
the stable operation of the cell [5]. Other factors that affect the 
temperature distribution and fuel utilization are the thickness 
and porosity of the anode. As a consequence, the geometrical 
design of SOFC is important in establishing a well-operating 
cell. In order to efficiently develop planar SOFC stacks, it is 
convenient and effective to have the capability to experiment 
numerically with the effects of different geometric designs on 
the operation and performance. 

In the past, modeling of the planar SOFC during steady 
operation to calculate the temperature and current density dis- 
tributions has been reported [6-12]. Investigations of planar 
SOFC operation and planar SOFC performance have predicted 
cell temperature and current distribution for various flow pat- 
terns [13-14]. However, only less work has been performed 
on the influences of variation of structure and operating con- 
ditions on the performance of SOFCs [15]. The purpose of 
present work is to demonstrate a model to predict tempera- 
ture distributions, species concentrations, current density and 
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Nomenclature 


Ck concentration of component k 
Cp specific heat (Jkg~! K7!) 


Dxett effect diffusivity of component k (m? s7!) 

Eç internal energy of mixture gas (J) 

Es internal energy of solid (J) 

F Faraday constant (C mol`!) 

I local current density (A/cm?) 

Ik mass source of mixture component k 

J transfer current density (A cm?) 

keff effect thermal conductivity coefficient (W m7!) 


kg thermal conductivity coefficient of gas (W m~ 1) 
ks thermal conductivity coefficient of solid (Wm!) 
K permeability coefficient (m?) 

P pressure of mixture (Pa) 

R gas constant (J mol`! K7!) 


Sk reaction coefficient of component k 
SE energy source (J) 

SM momentum source 

AS entropy change (J mol~!) 

T Temperature (K) 

U mixture velocity (m/s) 


Greek letters 


ô anode thickness (m) 

E porosity 

Nact.a activation potential at anode (V) 

Nact.c activation potential at cathode (V) 

Neonc concentration potential (V) 

Left effect viscosity 

p density (kg m~?) 

Oeff effect electrical conductivity (27! m7!) 


over potentials. In this study, co- and cross-flow cases were 
examined. In particular, the effects of the structural param- 
eters and the operating conditions on the performance of a 
SOFC are discussed in a co-flow case, such as thickness and 
porosity of anode, delivery rates and temperatures of fuel 
and air. 


2. Mathematics model 
2.1. Model geometry 


Generally, the repeating unit of a typical planar SOFC stack 
is constructed of a positive electrode—electrolyte-negative elec- 
trode (PEN) and an interconnector plate “stacked” together. For 
the sake of simplicity of calculation, one repeating cell unit, 
which includes a PEN, air and fuel channels and 1/2 intercon- 
nect thickness top and bottom in the Y-direction, was analyzed 
in this simulation. The one-cell stack and the single unit model 
are illustrated in Fig. 1. In this model, the thicknesses of anode, 
cathode, electrolyte and interconnect were 0.5, 0.25, 0.05 and 
1.0mm, respectively. 


DN 


1mm 


interconnector 


3mm 


Fig. 1. Illustrations of the one cell-stack and singe cell unit model. 


2.2. Thermo-fluid model 


The ANSYS-CFX code was selected to solve the thermo- 
fluid model. In the simulation, the solid and fluid domains were 
divided into some discrete meshes, and in each computational 
mesh, the conservation equations of species, mass, momentum 
and energy were solved using the finite volume method. 

In general, gas species transfer mainly by convection in the 
flow channels and diffusion in the porous electrodes. The species 
conservation equation: 


V(pCKU) = V(DxetfV Ck) + Ik, k = H2, O2, H20 (1) 


where J, is the rate of production or consumption of specie k, 
and given by [16]: 


Skj 
Ik = Ł+— 
N 


The diffusion coefficients of hydrogen and oxygen are 


obtained as follow [16]: 
nisl S. meal a 
Ta E 

(3) 


The mass conservation equation: 


(2) 


V -(eoU) = 0 (4) 


Both the air and fuel flows were considered as ideal gas 
mixtures with the density given by: 


_ P mk 
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where mm is the mass fraction of specie k with molecular weight 
My. 
The momentum conservation equation: 


2 ou 4 J p! ðw 
u w 
P ax dy Oz 
3u? x dv? P ðw? 
əx? ay 


Up Bice ( Ti © 
where Sm is momentum source, and SĮ =0 in the flow chan- 
nels. However, in the porous electrodes, Darcy law with constant 
porosity and permeability is applied to model the momentum 
source as follows [16]: 


Meff 2 
SM = —— E U 7 
M K (7) 


where eff is the effective viscosity of the mixture gas and is 
given by [16]: 


Xk Uk 
eff = = 8 
Leff FX (8) 
1 -y1/2 (M, / M014? 
RS [ + (uk/ 4j) ( 5 / k) ] (9) 


[8 + (M/M) 1/7 


where X, is the molar percent of the specie k, uj and upę are 
kinematical viscosities of specie j and k, respectively. 
The energy conservation equation 


V-(U(orEs + p)) + V (tU) + V - Kerr VT) + S&=0 (10) 


Heat transfer between the fluid and solid materials was limited 
to conduction and convection, and the effect of radiation was 
neglected in this calculation because it is very small relative 
to the other kinds of heat transfer. Additionally, the effective 
thermal conductivities of porous electrodes are calculated by 
the following Eq. (11) [17]: 


kefe = eke + (1 — £)ks (11) 


where kf and k; are thermal conductivities of fluid and solid, 
respectively. Sg is energy source term expressed by Eq. (12) 
and mainly consists of reaction and Ohmic heats [16]. 


SE = a 43 2an (12) 
E = Tk 5 IF Nact 


Specific heats of gas components are described as the functions 
of temperature [16]: 


Cp =atbx 10T +c x 10°T? (13) 


Table 2 
Properties of the solid material 


Table 1 

Coefficients of the specific heats of gas 

Gas a b c 
Hydrogen 25.8911 —0.8373 2.0138 
Oxide 29.0856 12.9874 —3.8644 
Water gas 30.3794 9.6212 1.1848 


where a, b and c are relevant coefficients, as listed in Table 1. 
Solid material properties used in this simulation are listed in 
Table 2. 


2.3. Boundary conditions 


The cell unit considered in this simulation represented a 
repeating unit in the middle of a large stack, and external walls 
of the cell unit were assumed to be adiabatic. Constant temper- 
ature, delivery rate, and gaseous composition were imposed at 
the inlet boundaries for the fuel and air. 


2.4. Electrochemical model 


2.4.1. Assumptions and reactions 
The oxidant reduction reaction occurring at the cathode is 
expressed as follows: 


102 + 2e7 > 0% (14) 


The oxygen ions transfer through the electrolyte and then into 
the active reaction areas of anode. The electrochemical reaction 
of fuel at the anode is: 


Hz + 07> > H20 + 2e7 (15) 
So the overall reaction is: 


Hz + 402 > H2O (16) 


2.4.2. Dynamics of electrochemical reactions 

According to Faraday’s law, the reaction rates depend on the 
current density i [7]: 
par aga (17) 

dt dt 

where df/dt, dO2/dt are the molar consumption rates of fuel and 
oxygen at the anode and the cathode, respectively. 

During the process of energy transforming, when the charge 
transfer reaction at the electrolyte—electrode interface is too slow 
to provide ions at the rate required by the demand of current, 


Cell component Density Effect thermal conductivity (W/m K) 
Interconnect 7700 13 

Anode 6200 6.23 

Cathode 6000 9.6 


Electrolyte 5560 2.7 


Specific heat (J/kg K) Porosity (%) Permeability coefficient (m?) 


0.8 - - 
0.65 25/35/45 1.0E—12 
0.9 35 1.0E—12 
0.3 - - 
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the activation polarization occurs and is defined by the Butler- 
Volmer equation [18]: 


2F 2F 
i = io {exp (-0 7a) = (ea = paras) } (18) 


Eq. (18) is simplified and described by the Tafel empirical 
formula [18,19]: 


RT . (19) 
o = l 
Nact,a 2Fioa 
RT . RT ; 
Nact,c = — (=) Inio,c + (=) Ini (20) 


where £ is the transmission coefficient and 6 =0.5 in this simu- 
lation, acta and Nact, are the activation potentials at the anode 
and the cathode, respectively. ig, and io, are the exchange cur- 
rent densities at the anode and the cathode, respectively. In this 
study, exchange current densities at the anode and the cathode 
are 5300 and 2300 A/cm?, respectively. 

The concentration polarization is also calculated as follows 
[18]: 


(21) 


Neon = 


RT n [zeae ad 
2F YHp,bulk YH20,A/E 


where YH, bulk and yH,0,bulk are the hydrogen and water 
gas molar fractions in the fuel flow channel, respectively. 
Ym, Aæ and yp, a are those in the interface of elec- 
trolyte/anode, respectively. 

Ohmic polarization is ignored in this simulation, so the open- 
circuit voltage is given by: 


V = E — Mact — Neon (22) 


where V is the open-circuit voltage and E is Nernst voltage 
obtained by: 


RT RT 1/2 
E = Ink + In (pra) x (Poa) 
(Pmo0) 


2 
2F 2F (23) 


where Ep is the standard voltage of the cell, pH,0, Po, and pp, 
are the partial pressures of water gas, oxygen and hydrogen, 
respectively. 


Table 3 


3. Simulation results and discussion 


In the calculations, the modeling tool couples an electro- 
chemical calculation method with a commercial computational 
fluid dynamics (CFD) simulation code. The finite-volume 
Navier-Stokes and transport equations are solved to obtain the 
gas species concentrations and temperatures at each position in 
the cell. The information is passed to the electrochemical model 
(subroutine). Then the local current density is calculated and 
applied to obtain the hydrogen reaction rate, heat source and 
species sources. Gas species concentrations and temperature 
distributions are then calculated for the next iteration, and so 
on, until convergence of solution is achieved. The cell operating 
conditions used for the simulation are listed in Table 3. 

Fig. 2 presents the temperature distributions of PEN for 
co-flow case and counter-flow case under the first working condi- 
tions. The average temperature of PEN is 979 °C with maximum 
and minimum temperatures of 1053 and 916 °C in co-flow case. 
Furthermore, it can be seen that the temperature of PEN increases 
uniformly along the direction of fuel flow, and is highest near 
the fuel outlet. However, for counter-flow case, average temper- 
ature is 996°C with maximum and minimum temperatures of 
1088 and 923 °C. In addition, it should be noted that the temper- 
ature of the PEN rises rapidly, reaching a maximum near the fuel 
inlet, and then gradually drops. Of the two flow configurations 
calculated, the co-flow case has the more uniform temperature 
distribution and a smaller temperature difference (A7) from the 
air inlet to the outlet. This is due to the offsetting effects of air 
near the inlet, at its coolest, being aligned with the fuel inlet. As a 
consequence, the temperature gradient is smaller in the co-flow 
case, which must result in a smaller thermal stress in the PEN 
although the thermal stress distribution was not investigated in 
this paper. 

Fig. 3 illustrates the hydrogen fraction distributions in the co- 
flow and counter-flow cases. Itis obvious that along the direction 
of fuel flow, the mass fraction of hydrogen on the interface of the 
electrolyte/anode decreases due to the electrochemical reaction. 
In particular, for the two flow case, the difference of hydrogen 
mass fraction is larger in co-flow case, so the utilization of fuel 
is higher accordingly. With increasing in inlet temperature of 


The cell operating conditions and parameters used for simulation 


Sample number Fuel 


Air Flow pattern 


Delivery rate Inlet temperature (K) 


Hydrogen mass 


Delivery rate Inlet temperature (K) 


(m/s) fraction (%) (m/s) 
1 0.5 973 0.8 3 873 Counter-flow Co-flow 
2 1.0 973 0.8 3 873 Counter-flow Co-flow 
0.5 
3 1.0 973 0.8 3 873 Co-flow 
1.5 0.8 
4 0.5 973 0.9 3 873 Co-flow 
1 3 
5 0.5 973 0.8 2 873 Co-flow 
1 
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Temperature [K] 
— 1.053e+003 


~= 1.038e+003 
M 1.023e+003 
1.007e+003 
9.923e+002 
9.77 le+002 
9.619e+002 
9.467e+002 
9.315e+002 


9.163e+002 
(a) 


Fuel flow 
Air flow 


Temperature [K] 
- 1.088e+003 
~ 1.070e+003 
m 1.051e+003 

1.033e+003 

1.015e+003 
9.968e+002 
9.785e+002 
9.603e+002 
9.421e+002 
9.239e+002 


Air flow 


Fuel flow 


(b) 


Fig. 2. PEN temperature distributions for the co-flow case (a) and counter-flow case (b) under first working conditions. 


Hydrogen Mass Fraction 
6.800e - 001 


6.156e - 001 
5.51 le - 001 
4.867e - 001 
4.222e - 001 
3.578e - 001 
2.933e - 001 
~ 2.289e - 001 
~ 1.644e - 001 
- 1.000e - 001 
(a) 


Fuel flow 
Air flow 


Hydrogen Mass Fraction 
6.800e - 001 


6.189e - 001 
5.578e - 001 
4.967e - 001 
4.356e - 001 
3.744e - 001 
3.133e - 001 
— 2.522e - 001 


-1.91le- 001 


Air flow 


Fuel flow 


-= 1.300e - 001 
(b) 


Fig. 3. Hydrogen mass fraction distributions on the interface of anode/electrolyte for the co-flow case (a) and counter-flow case (b) under first working conditions. 


the fuel gas, the average temperatures and temperature differ- 
ences (AT) also rise in the two flow cases, as shown in Fig. 4. 
Consequently, overall considering the temperature distribution 
and the hydrogen utilization, it is an advantage to choose the co- 
flow case for SOFC steady operation. For this flow case, other 
parameters that influence the performance of the SOFC are also 
studied. 


Temperature [K] 
= 1.099e+003 an 
~ 1.079e+003 

E 1.059e+003 
i 1.040e+003 
1.020e+003 
1.000e+003 
9.803e+002 
9.606e+002 
9.408e+002 


9.210e+002 
(a) 


I 


Fuel Flow 
Air Flow 


Firstly, we focus on the effects of the fuel gas on the temper- 
ature and current density. Fig. 5 shows the PEN temperature 
and the current density distributions. Increasing the delivery 
rate of fuel, the average temperature and temperature difference 
(AT) of PEN rise (see Fig. 5(a)), and the average current density 
also increases (see Fig. 5(b)). This indicates that it is effective 
for improvement of the electrical performance to increase the 


Temperature [K] 
Air Flow 


- 1.167e+003 
- 1.141e+003 

E 1.115e+003 
ler 1.089e+003 
1.063e+003 
1.038e+003 
1.012e+003 
9.858e+002 
9.599e+002 


9.340+002 
(b) 


Fuel Flow 


Fig. 4. PEN temperature distributions for the co-flow case (a) and counter-flow case (b) under second working conditions. 
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Fig. 5. Current density distributions on the interface of anode/electrolyte (a) and PEN temperature distributions (b) for the co-flow case under third working conditions. 
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Fig. 6. Current density distributions on the interface of anode/electrolyte (a) and PEN temperature distribution (b) for the co-flow case under forth working conditions. 


delivery rate of fuel. But the temperature gradient also rises, 
which may cause unexpected thermal stresses. Consistent with 
increasing of the hydrogen proportion from 80 to 90 and 100% in 
the fuel gas, the average current density, the temperature differ- 
ence and the average temperature also rise due to the increasing 
of reaction rate and more reaction heat being accumulated, as 
shown in Fig. 6. 

Figs. 7 and 8 indicate the effects of the fuel on the cell poten- 
tial. In the anode, the positive active potential directs the transfer 
current from the catalysts to the electrolyte, and vice versa in the 
cathode. It can be seen that the absolute value of the active poten- 
tial in the cathode is higher than that in the anode. Moreover, the 
active potentials at the anode and the cathode increase from the 
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fuel inlet to the outlet in the co-flow case, and the active potential 
is highest near the exit of the oxidant flow where it coincides 
with the regions of lowest oxygen concentration on the cathode. 
In particular, by increasing the fuel delivery rate and hydrogen 
mass fraction, the average active potentials at electrodes rise due 
to the air depletion. 

Then the influences of the air flow are also considered. In 
the SOFC system, air gas not only provides oxygen ions but 
also has a cooling function. Increasing the inlet velocity of air, 
the average and maximum temperatures of PEN drop, as shown 
in Fig. 9, which is because more reaction heat is absorbed and 
released by the air with the higher delivery rate, although the air 
utilization is reduced. 
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Fig. 7. Activation potential distribution at the anode (a) and cathode (b) for co-flow case under third kind of working conditions. 
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Fig. 8. Activation potential distributions at the anode (a) and cathode (b) for co-flow case under forth working conditions. 
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Fig. 9. PEN temperature distributions with air velocity 2 m/s (a) and 1.0 m/s (b) for the co-flow case under fifth working conditions. 


Finally, the effects of the anode structure on the concentra- 
tion potential of SOFC are considered. In this simulation, it is 
assumed that the electrochemical reactions occur at the inter- 
face of electrode and electrolyte. Since the diffusion of reactant 
from anode channel to the interface of anode/electrolyte is influ- 
enced by the thickness and porosity of anode, the concentration 
potential is also influenced, which must influence the output 
voltage. Fig. 10 illustrates the variation of the anode concentra- 


~ 
2 
= 


— a — anode thickness 0.25mm 
e — anode thickness 0.5mm 


— a — anode thickness 1.0mm 


Concentration potential (V) 


Distance from the gas inlet (dm) 


tion potential with the change of anode structure. It is observed 
that the concentration potential at the anode rises with the anode 
thickness. However, with the increase of anode porosity, the 
concentration potential drops. 

To validate the numerical simulation for the SOFC per- 
formances, some basic tests on the simulative model for the 
SOFC performances have been carried out, such as complex 
impedance testing (VMP/2, Princeton Applied Research Corp., 


o 
2 
o 
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Fig. 10. Concentration potential distributions at the anode with the variation of thickness (a) and porosity (b) of anode for co-flow case under second working 


conditions. 
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USA) and testing of the SOFC working performances (SF30, 
Wuhan Lixing Testing Devices Corp., China). Since the cal- 
culation conditions in this paper are different from the testing 
parameters, the present simulative results cannot be compared 
with the testing results. However, the present simulative results, 
such as the temperature distributions and the influence of the 
flow pattern on the temperature distributions are in qualitative 
agreement with the results in the reported literature [7,14]. 


4. Conclusions 


The exchange current density was used to couple the thermo- 
fluid model with the electrochemical model. Source terms were 
introduced in the governing equations to calculate the mass vari- 
ation of the reactants and products, and the momentum and 
energy in the reaction process. In particular, the Darcy law was 
applied to describe the flow in the porous electrode. The perfor- 
mances of the SOFC in the steady state were simulated and the 
conclusions were: 


(1) The co-flow case, due to its relatively uniform temperature 
distribution, smaller thermal gradients and higher average 
current density, offers thermo-structural advantages com- 
pared with the counter flow case. 

(2) For the co-flow case, with increase in the delivery rate of the 
fuel gas or hydrogen mass fraction in the fuel, temperature 
gradients in the PEN and the average current density rise, 
and the active potentials at the electrodes also rise, which 
may be harmful to the cell output power .... In an attempt 
to achieve a more uniform temperature distribution, it is 
effective to decrease the temperature gradients of the PEN 
by increasing the delivery rate of air. 

(3) Itis beneficial to decrease the concentration potential of the 
cell by decreasing the thickness and increasing the porosity 
of the anode according our simulation results. 
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